Velocity fluctuations of noisy reaction fronts propagating into a metastable state: 

testing theory in stochastic simulations 
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The position of a reaction front, propagating into a metastable state, fluctuates because of the 
shot noise of reactions and diffusion. A recent theory [B. Meerson, P.V. Sasorov, and Y. Kaplan, 
Phys. Rev. E 84, 011147 (2011)] gave a closed analytic expression for the front diffusion coefficient 
in the weak noise limit. Ifere we test this theory in stochastic simulations involving reacting and 
diffusing particles on a one-dimensional lattice. We also investigate a small noise-induced systematic 
shift of the front velocity compared to the prediction from the spatially continuous deterministic 
reaction-diffusion equation. 
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I. INTRODUCTION 



Effects of shot noise on the propagation of reaction 
fronts have been the subject of extensive studies in 
physics, chemistry and biology, see e.g. Refs. [l|, for 
reviews. The shot noise conies from the discreteness of 
reacting particles and random character of their reac- 
tions and migration in space. The role of shot noise in 
the front dynamics can be very different for fronts prop- 
agating into an unstable state and for those propagating 
into a metastable state (a state which is linearly sta- 
ble but non- linearly unstable). A front propagating into 
an unstable state can be extremely sensitive to the shot 
noise in the leading edge of the front [iHl]- In particu- 
lar, for pulled fronts, the front diffusion coefficient Df, 
describing random wandering of the front position, in- 
duced by the noise, is only logarithmically small with 
the typical number of particles A^ 3> 1 in the front re- 
gion |5|. Fronts propagating into a metastable state be- 
have more "normally" in this respect, and the scaling 
Df ^ 1/N was conjectured a while ago [JQ- Not much 
was known beyond this conjecture until recently, when 
a systematic theory of noisy reaction fronts propagat- 
ing into a metastable state was developed, and closed 
analytic expression for Df derived In this work we 
will test the theory of Ref. in stochastic simulations 
of the microscopic model of reacting and diffusing par- 
ticles on a one-dimensional lattice. In addition to the 
front wandering we will also address a more subtle effect 
of a systematic shift in the average front position, com- 
pared to prediction from the reaction-diffusion equation: 
a spatially continuum mean-field theory for this system 
@, Q ■ This systematic shift was previously addressed in 
Ref. , where it was suggested that it is exponentially 
small with respect to the carrying capacity of the system 
K, defined below. We revisit this issue here. 



II. THEORY 
A. General 

The microscopic model we will adopt here, see e.g. 
Refs. 0, H^, involves a single species of particles A re- 
siding on a one-dimensional lattice of sites. The number 
of particles on each site i varies in time as a result 
of two types of Markov processes [lH. The first type 
involves stochastic on-site reactions, such as A — >■ 2A, 
A ^ 0, etc. with given rates constants. Assuming for 
simplicity that, in each of these reactions, the number of 
particles changes by one, one can reduce all on-site reac- 
tions to a birth-death Markov process with given birth 
and death rates Xijii) and ii{ni), respectively. The sec- 
ond process is independent random walk of each particle 
between neighboring sites with rate constant Dq. 

A spatially continuous deterministic theory for this 
system can be derived under two assumptions Q. First, 
assuming that the carrying capacity is a large parameter, 
/•C ^ 1, one can write 



X{n,) = vKX{q.i) and /i(nj) = vK^i{q{), 



(1) 



where rii ^ 1, = rij/ K, X{qi) ~ fiili) ~ 1, and ly 
is a rate constant [l^iliSl- When neglecting noise, this 
brings about a spatially-discrete version of the determin- 
istic reaction-diffusion equation [l^ 



<}i = ^fili) + Da{qi-i + qi+i - 2qi) 



(2) 



where f{qi) = X{qi) — fi{qi) is the rescaled birth-death 
rate function. Assuming, in addition, a sufficiently fast 
hopping, Dq 3> ly, one can use a continuous position 
coordinate x instead of the discrete index i and transform 
Eq. ^ to the reaction-diffusion equation d, Q 



dtq^uf{q) + Ddlq, 



(3) 



where D = DqH^ is the diffusion coefficient of the parti- 
cles, and h is the lattice constant. 

A spatially continuous deterministic reaction front cor- 
responds to a traveling wave solution (TWS) of Eq. ([3]), 



2 



q{x,t) = qoiOy ^ — ^ ~ ^0^5 that connects two linearly 
stable homogeneous states described by zeros of the func- 
tion /(g). Let one of them, g = > 0, be located at 
X — >■ — oo, and the other, g = 0, at a; — > oo. There is also 
a linearly unstable state g^, so that < g„ < g*. The 
TWS qo{^) obeys the equation 



Dq'o+c„q'„ + l^f{q„)^0 



(4) 



and is unique up to an arbitrary shift in ^ (1, Q . Let us in- 
troduce the effective potential V{q) = f{u) du. Then 
the state g = is retreating (co > 0) when ^(g*) > 
d, Q . In this case g = is called the metastable state. 
For V(g*) < the state g = g* is retreating (cq < 0), 
and so it is metastable. Finally, for ^(g*) = the de- 
terministic front is standing, cq = 0. We will call the 
corresponding set of parameters the stall point. The typ- 
ical front width is on the order of the diffusion length 
Id = {D/v)y\ 

To account for shot noise, one can go back to the 
stochastic model on the lattice and use the master equa- 
tion describing the evolution of the multivariate proba- 
bility distribution P(n, t) = P(ni, n2, . . . ,t). For typical, 
small fluctuations this master equation can be approxi- 
mated, via a truncated expansion, by a Fokker-Planck 
equation 0, Instead of dealing with the Fokker- 

Planck equation directly, Meerson et al. derived a 
Langevin equation to which the Fokker-Planck equation 
is equivalent. For Dq ^ this Langevin equation is 
continuous in space and takes the form 0] 



dtq{x,t) = yf{q) + Ddlq^ 
+ 9. 



ug{q) h 



2q{x,t) Dh 



K 



K 



ri{x,t) 



(5) 



where g{q) — A(g) -I- p,{q) is rescaled on-site diffusion 
coefficient in the space of population sizes, whereas ri{x, t) 
and x(x, t) are independent Gaussian noises which have 
zero means and are delta-correlated in x and in t: 

t)7i[x',t')) = {x[x, t)x{x',t')) = 5{x - x') S{t - t'). 

(6) 

The (multiplicative) noise terms in Eq. ([S]) come from the 
on-site birth and death processes, and from the random 
walk (the third and fourth terms on the right, respec- 
tively) . 

Rescaling the position coordinate x and time t hy Id 
and l/iy, respectively, one can rewrite Eq. ([S]) as 

dtqix, t) = f{q) + dlq{x, t) + N-^'^R{x, t, q) , (7) 

where 

R{x,t,q) - ^^i^{x,t)+d^[^qx{x,t)]. (8) 

As the parameter e = N^^^^ ^ 1 is small, one can solve 
Eq. ([7]) perturbatively around the deterministic TWS 
with an a priori unknown front position, slowly varying 



in time. The first order of thisperturbation theory yields 
a closed-form analytic result Q for the front diffusion co- 
efficient Df in the frame moving with the average front 
velocity c. In the dimensional variables 



Df = 



D 



(9) 



where 



So 



I-oo ('Z^e^o?)2 5(go) + go [(g(,e^o?)'] } 



(10) 

It follows from Eq. ([9]) and the definition of N that Df 
scales as with the carrying capacity parameter K, 
and as _D^/^ with the particle diffusion coefficient D. 

The probability distribution of the typical position X 
of the fluctuating front at time r is a Gaussian: 



P(X,r) 



\ATr Dt 



1/2 



exp 



soN{X - erf 



ADt 



(11) 



The Gaussian asymptote holds in the region of \X/t — 
c\ < (vD)^/^ H. Predictions ([9l)-(fTT|) were also ob- 
tained from a WKB approximation to the master equa- 
tion, in conjunction with linearization of the WKB equa- 
tions around the deterministic TWS solution 0. 

In the first order of the perturbation theory in e, that 
the authors of Ref. Q confined themselves to, the aver- 
age front velocity does not change compared with the pre- 
diction of the spatially continuous deterministic theory: 
c = Cq. In the second order in e = N~^/'^ <C 1 a small 
systematic shift 5c of the average front velocity should 
appear, caused by the shot noise: c = cq + 5c. In the 
rescaled variables 5c is expected to behave as (Sc = xe^, 
where x is a dimensionless factor, which magnitude and 
sign are determined by the properties of functions /(g) 
and 5(g), see Eqs. ^ and ([8]). In dimensional variables 



5c = xf-Dve = ——. 



(12) 



That is, the shot-noise-induced velocity shift is expected 
to scale as . This is in contrast to the scenario of Ref. 
[lo| which predicts a front velocity shift exponentially 
small in K. Somewhat surprisingly, the front velocity 
shift (|12p is independent of D. The actual second-order 
perturbative calculations are very cumbersome, and we 
will not embark on this route here, confining ourselves to 
the simple argument leading to Eq. (jl2[) with the unde- 
termined dimensionless factor x. 

What is the expected validity range of Eqs. ([9]) and 
([T^? Equation ^ only requires I?o ^ v (fast hopping) 
and K ^ 1 (weak shot noise). Equation (|12p demands 
an additional criterion. This is because there is an addi- 
tional, well-known mechanism of shift of the average front 
velocity. It is purely deterministic and comes from the 
discreteness of the lattice, see Eq. ^ [l3|. At small v/Dq 
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this velocity shift behaves as {Sc)d = const x {v/Dq)cq, 
as shown by Keener This deterministic velocity 

shift will mask the shot-noise-induced velocity shift un- 
less \5c\ ^ \{Sc)d\- This criterion can be rewritten as the 
following strong inequality: 



Doh 
Kcq 



> 1. 



(13) 



Far from the stall point the front velocity cq can be 
roughly estimated as cq ^ i^Id, and the two criteria, 
K ^ 1 and Eq. become a strong double inequal- 
ity 



1 « K « {Do/i^y/^ 



(14) 



that is satisfied for sufficiently fast hopping, Do/iy ^ 1. 
Close to the stall point cq is small, and the criterion 
is easily satisfied, making the right inequality in Eq. 
Plj) unnecessary (one still needs to demand Dq » v). 
This feature makes standing fronts, cq = 0, especially 
convenient for numerical studies of small front velocity 
fluctuations and shift induced by shot noise. Another 
(even more apparent) advantage of the case of cq = 
is the possibility to perform the simulations in relatively 
short systems. 



B. Example 

A simple set of on-site reactions which exhibits bista- 
bility is A — >■ and 2A 3A. Let the rate constants 
of these three processes be fiQ, Aq and ctq, respectively. 
Following Ref. 0, we will perform rescaling so that 
X{q) = 2q^ and Ji{q) = 717 + 5^, where v ~ 3Aq/(8(To), 
K = 3Ao/(2(Jo), and 7 = 8/iocro/(3AQ). The on-site dy- 
namics exhibits bistability at (5^ = 1 — 7 > 0. In this case 
the zeros of function f{q) = —q{q — qu){q — q*) describe 
two stable fixed points of the deterministic theory, and 

= 1 + (5, and an unstable fixed point q^ = 1 — 5 such 
that < < q*. Here the TWS <Zo(0 elementary 



90(e) 



Co 




According to the spatially continuous deterministic the- 
ory, the state g = is advancing at 1/3 < 5 < 1, re- 
treating atO < 6 < 1/3 and stationary at S — 1/3. 

The function g{q), entering Eq. (|10p . is equal to g{q) = 
A(q) + M(g) =7g-f 2g2 + g3_ 

In this example, the factor sq = so(6) from Eq. ^TU\\ 
can be evaluated analytically for any (5 [7| . In its turn, the 
factor X from Eq. (fT2|) also depends only on S, but this 
dependence is unknown at present. In a similar bistable 
system, studied in Ref. [Ifll, the front velocity shift Sc 
even changed sign when a parameter similar to our 6 was 
varied. In the present work we focus on the particular 
value of S = 1/3, corresponding to the stall point of the 



mean field theory Here sq = ^2/6 0, Ao = 9fio/{2K), 
ctq — 27 fiQ / {AK'^) and 1/ = 9^q/8, and we can rewrite 
Eqs. (HSI, ®, (HI]) and (IT2|) as 



Df 
P(X,r) 

5c 



3{l + eV 

2K ' 
1 

y^AnDfT 

"IT' 



^^x~Sct, (16) 



exp 



ADfT 



(17) 
(18) 

(19) 



where Eq. (|19p includes an unknown numerical coefficient 
a = 0(1). It is Eqs. (|16p - (fT9|) that we aimed at verifying 
in stochastic simulations. 



III. SIMULATIONS 

We performed extensive Monte Carlo simulations on a 
one-dimensional lattice with spacing h — 1 and integer 
length L ^ 1. Every lattice site < z < L can be 
occupied by any number of particles n^. The number of 
particles changes, at each step, because of independent 
random walk or one of the reactions A — >■ and 2 A ^ 3^, 
see Sec. Ill Bl We employed the Gillespie algorithm [l6l |. 
A site J is chosen with probability proportional to rij. 
Then a particle on that site is chosen at random. It can 
jump to a neighboring site on the right or left, produce 
a new particle, or die with probabilities 

Pright = Pleft = D / {bl + di + 2D), (20) 
Pbirth = 61/(61+^1+21?), 

Pdeath = di/{bi+di+2D), 
respectively. Here 

61 = Ao(7ij — l)/2 and rfi = /.io + croinj — l)(?^j — 2)/6. 

After every single-particle event, the time is advanced 
by l/[M(6i +di + 2D)], where M = J^j is the total 
number of particles in the system. At the boundaries 
i = and L, particles are not allowed to migrate out 
of the system, which corresponds to no-flux boundary 
conditions. We took the initial number density profile 
that corresponds to the mean-field front solution given 
by Eq. (fTB]) : nj{t = 0) = K qo{j), such that the front 
interface is located in the center of the system, Xq = L/2. 

We performed many simulations and determined the 
front position for different values of parameters. The 
number of simulations for each value of K and D varied 
between 60 and 420. A typical density profile Uj, ob- 
served in the simulations, is shown in Fig. [T] To find the 
front position X, we fitted the density rij by the mean- 
field profile that moved distance X — Xq; the value of 
X served as the fitting parameter. For a fixed simula- 
tion time T, stochastic fronts corresponding to different 
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FIG. 1: (Color online). A typical realization of the front den- 
sity profile observed in a Monte Carlo simulation of the lattice 
model with reactions A — ^ and 2A ^ 3A and unbiased ran- 
dom walk. The smooth line shows the mean-field solution 
n{j) = K qo{j), see Eq. (|16|l . with a fitted front position. The 
parameters are: K = 30, D = 20, v = 0.1, S = 1/3, and 
ft = 1. 



realizations propagate a different distance X — Xq. We 
computed the standard deviation ast of the distribution 
of the front positions X, and determined the fronts dif- 
fusion coefficient Df = a1^/{2T). Figures [5] a and b show 
the simulation results for Df as a function of the carrying 
capacity K and of the diffusion coefficient of the particles 
D, respectively. These plots also show theoretical predic- 
tions from Eq. (jl7p . without any fitting parameters. The 
agreement is good in both cases. 

Furthermore, we combined our sets of measurements 
for different values of parameters to obtain a rescaled 
probability distribution of the front positions, and there- 
fore of the average front velocities. First, in each set of 
simulations with the same parameters, we transformed 
to a moving frame, X' ~ X — X ^ hy shifting the mea- 
sured distribution by the average displacement X of the 
fronts in this set of simulations. Then we combined the 
data from the different sets of simulations on a single 
plot that shows the rescaled probability distribution of 
X', P{y) — yjADfTP{X'), versus the rescaled front po- 
sition 



y 



X' 



sqN 
ADt 



X' 



see Eq. (|T8| . This plot is shown in Fig. [3l As one can 
see, there is a very good collapse of data for different sets 
of simulations and a good agreement with the theory. 

We have also measured the average shift in the front 
position X, with the aim of verifying the scaling rela- 
tion (|19p for Sc. To remind the reader, for S = 1/3, the 
front velocity, predicted by the mean field theory, is zero. 
The measured average shift, therefore, comes from the 
shot noise. Figure |4] shows the front velocity correction 
Sc = X/t as a function of the particle diffusion coeffi- 
cient D for three different carrying capacities K. One 
can see that, at sufficiently large D, Sc becomes indepen- 
dent of _D, as predicted by Eq. This equation also 
predicts that, in the large-Z? limit, the velocity corrcc- 
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FIG. 2: (Color online). The front diffusion coefficient Df as 
a function of K (a) and D (b). Symbols with error bars: 
the simulation results. The solid lines are computed from 
Eq. ^7}. Theparamet ers are: D = 40 (a) and K = 30 (panel 
b). Also, u — 0.1, S = 1/3, the simulation time is r = 2000, 
and the system size is L = 600. 



0.6 
0.5 
0.4 
J 0.3 
0.2 
0.1 





FIG. 3: (Color online). Symbols: the rescaled distribution 
of front positions measured in simulations for different values 
of D and K (the total of 1591 simulations). The solid line: 
the normal distribution P{y) = tt"^/'^ exp(— y^) predicted by 
theory 0] for typical fiuctuations. No fitting parameters are 
used. The rest of parameters are: v = 0.1, 5 = 1/3, the 
simulation time is r = 2000, the system size is L = 600. 



tion must decrease with K as 1/K. This prediction is 
fully supported by our measurements, see Fig. [SJ where 
the rescaled quantity KSc is plotted versus D, and good 
collapse of data points for K = 60, 90 and 120, and a 
single data point for K = 150, is observed. For the case 
oi S = 1/3, considered here, the coefficient a in Eq. (E 
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is about 7.6. 
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FIG. 4: (Color online). The noise-induced velocity shift 
5c = X /t as a. function of D for three values of the carrying 
capacity: if = 60 (circles), if = 90 (asterisks) and K = 120 
(squares). The rest of parameters are the same as in Fig. [S] 
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IV. SUMMARY 



Our stochastic simulations strongly support theoreti- 
cal prediction of Ref. of the diffusion coefficient Df 
of noisy reaction fronts propagating into a mctastable 
state. The inverse dependence of D/ on K and the ^/D 
dependence oi Df on D are both verified, as well as the 
Gaussian distribution of the (typical) front positions in 
the reference frame moving with the average front speed. 

The simulations also clearly show that, for sufficiently 
high hopping rates, the noise-induced systematic shift of 
the front velocity is a second order effect in the small 
parameter e = 1/ViV, and so it scales as K~^. 

The front diffusion coefficient comes from the Gaussian 
region of the fluctuations of the front position. These 
are typical, small fluctuations. The WKB theory of Ref. 
also dealt with rare large fluctuations and predicted 
non-Gaussian tails of the front velocity distribution. Ob- 
serving these tails in stochastic simulations requires a 
more specialized simulation algorithm and is left for fu- 
ture work. 
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FIG. 5: (Color online). The rescaled noise-induced velocity 
shift KSc as a function of D for four values of the carrying 
capacity: K = 60 (circles), 90 (asterisks), 120 (squares) and 
150 (a single point, denoted by diamond). The rest of param- 
eters are the same as in Fig. O The dashed line corresponds 
to KSc = 0.68. Then, using Eq. (EH) with /io = (8/9)i^, we 
obtain a ~ 7.6. 
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